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Abstract 

A finite element approach for integrated fluid- 
thermal-structural analysis of aerodynamically heated 
leading edges is presented. The Navier-Stokes 
equations for high speed compressible flow, the 
energy equation, and the quasi-static equilibrium 
equations for the leading edge are solved using a 
single finite element approach in one integrated, 
vectorized computer program called LIFTS. The fluid- 
thermal-structural coupling is studied for Mach 6.47 
flow over a 3-inch diameter cylinder for which the flow 
behavior and the aerothermal loads are calibrated by 
experimental data. Issues of the thermal-structural 
response are studied for hydrogen cooled, super 
thermal conducting leading edges subjected to intense 
aerodynamic heating. 

Nomenclature 

A coolant passage area, eq. (5), or finite 

element area, eq. (11) 

c specific heat, eqs. (4) and (5) or fictitious 

damping constant, eq. (7) 

E,F x and y flux components 

H load vector, eq. (5) 

h convective heat transfer coefficient 

k thermal conductivity 

I, m components of unit normal vector 

m coolant mass flow rate 

[M] mass matrix 

[N] element interpolation function 
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coolant passage perimeter 

pressure 

heat flux 

load vector 

temperature 

reference temperature for zero stress 
surrounding medium temperature 
time 

time step 

conservation variable 

flow velocity components, eq. (2), or 
displacement components, eq. (7) 

coordinate directions 

density 

flow total energy, eq. (2), or strain 
components, eq. (18), or emissivity, eq. 

( 19 ) 


o Stefan-Boltzmann constant 

o x , Oy, txy fluid stress components, eq. (2), or solid 
stress components, eq. (7) 


Subscripts 
F fluid 

T thermal 

s structural 

Superscript 

n time step index 
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Introduction 

Design of lightweight structures and thermal 
protection systems for hypersonic cruise and reentry 
vehicles depends on accurate prediction of the 
aerothermal loads, structural temperatures and their 
gradients, and the structural deformations and 
stresses. Traditionally, an aerodynamicist will predict 
the surface pressures and heating rates assuming a 
rigid isothermal body. These aerodynamic heating 
rates are used by a structural heat transfer analyst to 
predict the structural temperature distribution. Finally, 
a structural analyst uses the temperature distribution 
and aerodynamic pressures to predict the structural 
deformations and stresses. Such traditional 
independent approaches require several iterations 
between the different analysis methods and analysts. 
The approach is relatively inefficient because the 
incompatible mathematical models require extensive 
postprocessing to transfer data. Moreover, the 
interdisciplinary coupling and interactions, which are 
not insignificant, are rarely attempted because the 
iterative process not only requires several additional 
solutions but also remodeling in each analysis. The 
coupling occurs primarily through the thermal 
response of the structure because: (1) the surface 
temperature affects the external flow by changing the 
amount of energy absorbed by the structure, and (2) 
the temperature gradients in the structure result in 
structural deformations which alter the aerodynamic 
surface and hence the flow field and attendant surface 
pressures and heating rates. Hence, an integrated 
interdisciplinary analysis procedure that would provide 
accurate, timely prediction of the coupled response is 
highly desirable. 

The interdisciplinary coupling and interaction has 
been demonstrated in the design of high speed 
vehicle structures. One such structure, a metallic 
Thermal Protection System (TPS) panel 1 for reentry 
vehicles, thermally bows into the airstream to relieve 
thermal stresses. The bowing creates an aerodynamic 
surface consisting of arrays of spherical like 
protuberances. As a first step to understanding the 
effect of the structural deformation on a hypersonic 
flow field caused by the interactions, the aerothermal 
loads on a dome surface in a Mach 6.5 stream were 
predicted computationally 2 and measured 
expenmentally 3 . These studies showed that the dome 
configuration caused complex shock and expansion 
waves and attendant localized heating rates that were 
2 to 6 times higher than those for a uniform flat surface. 
The existance of the high localized heating rates 
indicates that several more iterations between 
interdisciplinary analyses would be required to 
converge the aerothermostructural response of the 
TPS. Typically, these additional calculations were 
never performed. 

To further understand the coupled fluid-thermal- 
structural interactions, a thin panel representing a 
typical sidewall of an actively cooled engine structure 
has been analyzed with an integrated fluid-thermal- 
structural procedure 4 - 5 . The coupled analysis showed 
that the aerodynamic flow field was altered: (1) as the 
initial uniform wall temperature became nonlinear, and 


(2) as the panel deformed. The nonlinear temperature 
distribution and panel deformation not only changed 
the stress distribution in the panel, but also induced a 
complex external flow phenomena consisting of 
separated flow with shocks and flow expansions. The 
study confirmed that the interactions of the flow 
behavior, the panel temperature, and the panel 
deformation should not be neglected and that the 
interactions between the three disciplines could be 
coupled into a single analysis method. 

Leading 1 edges for hypersonic vehicles that 
experience intense stagnation point pressures and 
heating rates are often a challenge to the designer. 
For engine leading edges, such as cowls or the fuel 
injection strut shown in Fig. 1 , the intense loads can be 
amplified by an order of magnitude 6 when the leading 
edge bow shock is impinged upon by an oblique 
shock wave. The intense localized heating causes 
severe temperature levels and gradients. This makes 
the analytical prediction difficult because of the strong 
nonlinearity that is present in all phases of the thermal 
and structural analyses. The strong nonlinearity raises 
uncertainity in applying the existing analysis approach 
to such highly nonlinear thermal-structural problems. 
However, nonlinear behavior is important in the design 
of a viable structure under such severe conditions, 
hence the requirements of such an analysis procedure 
needs to be better defined. 

The purposes of this paper are to: (1) study the 
coupled fluid-thermal-structural behavior for 
aerodynamically heated leading edges using the 
technique developed for the heated panels, and (2) 
identify the critical issues or requirements for 
integrated fluid-thermai-structural analysis. The 
solution of the Navier-Stokes equations for predicting 
aerodynamic heating and the solution for the 
associated thermal-structural equations are obtained 
using a Taylor-Galerkin algorithm in one integrated, 
vectorized computer program called LIFTS (the 
Langley Integrated Eluid Ihermal Structural analyzer). 

The fluid-thermal-structural formulation used in 
LIFTS and the solution approach are described. The 
integrated analysis will be demonstrated on a three 
inch diameter cylinder used to obtain the experimental 
pressures and heat transfer rates to which the 
predictions will be compared. The issues associated 
with nonlinear thermal-structural response due to 
intense aerodynamic heating will be demonstrated by 
studies of a hydrogen cooled leading edge subjected 
to shock wave interference heating. The latter 
example will include the issues associated with the 
effects of severe temperature gradients on material 
property variations and thermal stresses. 

Flukf-Thermal-Stnj ctural Formulation 


The equations for aerodynamic flow are described 
by the conservation of mass, momentum, and energy 
equations. These equations can be written in the 
conservation form 


Aerodynamic Flow 
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where the subscript F denotes the flow analysis, {U F } is 
the vector of the conservation variables; (E F ) and {F f } 
are vectors of the flux components in the x and y 
directions. These vectors are given by 
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where p is the fluid density, u and v are velocity 
components, and e is the total energy. Each of the flux 
vectors contains two vectors of flux components 
representing the inviscid and viscous flux components. 
In the inviscid flux components, the pressure p is 
related to the total energy assuming a calorically 
perfect gas (constant ratio of specific heats). In the 
viscous flux components, the stresses <r x , o y , and x xy 
are related to the velocity gradients assuming Stokes' 
hypothesis 7 . The heat fluxes q x and q y are related to 
the temperature gradients by Fourier's law. The 
temperature-dependent viscosity is computed from 
Sutherland's law 7 and the thermal conductivity is 
computed assuming a constant Prandtl number of 
0.72. 

Structural Heat Transfer 


coolant temperature, (2) a convective heat transfer 
coefficient that represents the thermal conductance 
between the structure and the coolant, and (3) the 
coolant mass flow rate 8 . The energy equation for the 
coolant flow in the local x direction can be written in 
the form of eq. (3) where 

U T = p ft T f 

E t = mCjTj/Aj - kjdTj/dx 

F t = 0 (5) 

Hr = hP(T -T f ) 

In this equation, p f is the coolant density, c< is the 
coolant specific heat, T f is the coolant bulk 
temperature, m is the coolant mass flow rate, and kf is 
the coolant thermal conductivity. Therefore, Ej 
consists of the energy transport by convection (1st 
term) and conduction (2nd term), and Hx represents 
the heat transfer between the structure and the 
coolant. 

Structural Response 

The structural response is described by the quasi- 
static equations of motion which can be written in 
conservation form 

f (U s ) + f (E ) + f (F ) = 0 (6) 

dt dx dy 

where the subscript s denotes the structural analysis, 
{U s } is the displacement vector; {E s } and (F s ) are 
vectors of the stress components and are given by 


The thermal response of the structure is described 
by the energy equation which can be written in 
conservation form 

i (Ut)+ t (Er)+ 1 (Ft)=Ht (3) 

where the subscript T denotes the thermal analysis, Uj 
is the conservation variable, E T and F T are the flux 
components, and Hx is the heat load. For transient 
heat conduction without internal heat source, these 
quantities are 

Ut= p s c s T s 
Ex^dx 
F T = q y 


(U 8 ) T = [cu 

cv ] 



V 

(7) 

T 

(F } =[-T 
1 i J *y 

°y ] 



where u and v are displacement components in the x 
and y directions respectively; c is a psuedo-damping 
constant used to facilitate marching to a steady-state 
quasi-static solution. The panel stress components a x , 
<r y , and x xy are related to the displacement gradients 
and the temperature by the generalized Hooke’s law. 

Solution Procedure 


H T =* 0 (4) 

where p 8 is the material density, c s is the material 
specific heat, and T s is the temperature. The heat 
fluxes q x and q y are related to temperature gradients 
by Fourier's law. 


An explicit time marching finite element scheme, 
the Taylor-Galerkin algorithm, described in Refs. [4-5] 
is used to solve the fluid thermal structural equations 
(1) - (7). For brevity, only the essential features of the 
algorithm are highlighted herein. 

Taylor-Galerkin Algorithm 


For a convectively cooled structure, the thermal 
mass transport by the coolant is prescribed by an 
energy conservation equation based on: (1) a bulk 
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The basic concept of the Taylor-Galerkin algorithm 
is to use: (1) Taylor series expansion in time to 



establish recurrence relations for time marching, and 
(2) the method of weighted residuals with Galerkin's 
criterion for spatial discretization. The fluid, thermal, 
and structural conservation equations are written in the 
form of a scalar equation, 


dU dE 3F „ 


( 8 ) 


The key feature of the algorithm is to express the 
variation of the element fluxes E and F in the same 
form as the element dependent variable U, that is 


U(x.y.t) = [N(x,y)]{U(t)) 

E(x,y,t) = (N(x,y)](E(t)] ( 9 ) 

F(x.y.t) = {N(x.y)]{F(t)} 


where [N(x,y)J denotes the element interpolation 
functions, and {U}, {E}, and {F} are the vectors of the 
element nodal quantities. 

For flow analysis, the computation proceeds 
through two time levels, t n+1 /2 and t n+ i. At time level 
t n+ i/ 2 . values of U are assumed constant within each 
element. At time level t n+ i, the nodal values of U are 
computed assuming a linear element variation. The 
final equations are in the form, 

(Ml {AU F ) n+ ' = {Rp}" + (RpJj + {Rp}^ (10 > 


where AU n+1 => U n+1 - U n and [M] denotes the mass 
matrix, 


[M|=J{N}[NldA {11) 

A 

The first two vectors on the right hand side of eq. (10) 
represent the viscous flux over the element area and 
along the element boundary, 


IRp), = -At ( 


fi— 

i 3 x 


(N]dA (E) r 



— )[N]dA{F) n ) 
dy 


( 12 ) 


For thermal analysis, a single time level version of 
the algorithm is used 10 . Derivation of the element 
equations proceeds in the same fashion as described 
for the flow analysis. The final element equations are 
in the form, 

[M] (AU T ) n+1 = (R^* (R t )^ ( 14 > 

where 


(R T )"=At(| (_ }[N]dA{E r ) n 

a dx 
A 

+ f { — } [N] dA (F T ) n ) 05 ) 

A ^ 


(R^ = -At J(N) [N] dS (1 (E T )" + m (F T ) n ) 
s 

- -At J (N) [N] dS (q) " ( 16 ) 

s 

The vectors (Ej) and (F T ) contain the element nodal 
heat flux components. The vector (q) in eq. (16) 
represents nodal heat fluxes normal to the element 
surface boundary. Typical nodal heat flux 
components, for an isotropic material as an example, 
are computed from Fourier’s law, 


E^q^-kCD- 


F T =q y = -k(T) 


£T 

dx 

3 T 


( 17 ) 


3 y 

The nodal temperature gradients 3T/3x and 3T/9y 
depend on element types and element nodal 
temperatures, and can be computed explicitly. 
Benefits of the Taylor-Galerkin algorithm include: (1) 
all element matrices and element nodal quantities can 
be evaluated in closed form, i.e. , numerical integration 
is not required, even for quadrilateral elements, (2) 
material nonlinearity such as temperature dependent 
thermal conductivity can be included directly (eq. (17)), 
(3) nonlinear boundary conditions which will be 
described in the subsequent section can be 
implemented easily, and (4) an explicit time marching 
scheme can be used for transient nonlinear analysis, if 
preferred, thus avoiding the iterative solution of a set 
of simultaneous equations. 


(Rp)j = AtJ{N) |N) ds ( 1 ( E) n + m{ F) n ) (13) 
s 

where At is the timestep and the coefficients I and m 
are the components of a unit vector normal to the 
boundary. The last vector on the right hand side of eq. 
(10) represents both the inviscid and viscous fluxes 
evaluated at time level t n+ i/ 2 - All matrices in eqs. (11) 
(13) are evaluated in closed form 9 . To produce an 
explicit algorithm, the mass matrix [M] in eq. (11) is 
diagonalized. 


The same approach used for the thermal analysis 
is applied to derive the finite element equations for the 
structural analysis 10 . The equations are identical to 
eqs. (14) to (16) except that the subscript T is 
replaced by the subscript s everywhere. The vectors 
(E s ) and (F s ) now represent the element nodal stress 
components. For an elastic orthotropic material, 
typical nodal stress components in two dimensions are 
obtained using constitutive relations, 

a-c^-^Cr-T) ij = 1 . 2,3 ( 18 ) 
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where T Cl is the reference temperature for zero thermal 
stress. The material elastic constants cy and the 
thermal expansion parameters ft may be temperature 
dependent. Large-strain displacement relations are 
available in LIFTS and can be included in the 
computation of the strain components ej. 

The explicit Taylor-Galerkin procedure described 
herein is conditionally stable. The allowable time 
steps depend on the size of elements and the 
characteristics of the problem. Guidelines for 
determining allowable time steps have been 
established and can be found in Refs. 10-11. 

Initial and Boundary Conditions 

The fluid, thermal, and structural eqs. (1), (3), and 
(6) are solved subject to appropriate initial and 
boundary conditions. The initial conditions for these 
problems consist of specifying the distributions of the 
conservation variables {U} at time zero. 

The boundary conditions for supersonic flow 
consist of specifying all conservation variables along 
the in-flow surfaces. On supersonic outflow surfaces, 
the finite element formulation provides appropriate 
boundary conditions. A no slip condition (velocities set 
to zero) is specified at the solid surface. 

The boundary conditions for the thermal analysis 
are applied via the vector shown in eq. (14) and 
expressed in eq. (16). The surface nodal heat flux q is 
replaced by the quantities representing different types 
of thermal boundary conditions, 


r o 


(insulated) 


<1 = 


< 


h (T-TJ 


(specified heating) 
(surface convection) 


ea(T^-T^) 


(surface radiation) 


(19) 


The boundary conditions for the structural analysis, 
such as the applied pressure, can be added into the 
structural equations via the surface boundary vector. 
The procedure is identical to that for the thermal 
analysis previously described and is therefore omitted. 

For fully integrated analysis, the aerodynamic 
pressure and surface heating rate are integral parts of 
the analysis and are not boundary conditions. The 
Taylor-Galerkin approach permits the different 
disciplines to be combined easily for the coupled 
analysis. As an example, the coupled solution for the 
flow and the structure temperature can be obtained by 
solving the equation that combines the flow energy eq. 
(10) and the structure energy eq. (14). The result 
contains an interface temperature distribution valid for 
both the flow field and the structure and provides the 
continuity of heat flux across the boundary. 


Applications 

Two applications are presented to study the fluia- 
thermal-structural behavior of aerodynamically heated 
leading edges. A three inch diameter stainless steel 
cylinder subjected to aerodynamic heating from a 
Mach 6.47 flow is used as the first example to 
demonstrate the integrated analysis and validate 
predicted surface pressure and heating rates with 
experiment. The second example is a study of the 
thermal-structural response of a 0.25 inch diameter 
hydrogen cooled leading edge subjected to intense 
shock wave interference heating. Critical issues for 
integrated fluid-thermal-structural analysis and design 
of the leading edges are identified. 

Mach 6.47 Row Over a Cylinder 

The solution to Mach 6.47 flow over a cylinder is 
used to demonstrate the integrated flow-thermal- 
structural analysis approach and validate the 
analytical prediction of the aerothermal loads. A 
schematic of the experiment, performed in the NASA 
Langley 8-foot High Temperature Tunnel, is shown in 
Fig. 2. A 3-inch diameter, 0.5 inch thick, stainless steel 
cylinder was mounted on the panel holder and 
subjected to a uniform high enthalpy Mach 6.47 flow 
Approximately 50 pressure taps and coaxial 
thermocouples were placed circumferentially along the 
cylinder surface to accurately determine the 
aerodynamic pressure and heating rate distribution. 
Details of the experimental configurations, the tunnel 
flow conditions, and the experimental results are given 
in Ref. 6. 

The finite element model representing the flow 
domain and the cylinder is shown in Fig. 3. Due to 
symmetry, only one-half of the incoming flow domain 
and one-fourth of the cylinder are modeled. The flow 
field is characterized by the bow shock that stands off 
the cylinder and the thin boundary layer at the cylinder 
surface. Sharp gradients in the flow variables occur in 
these regions, hence closely spaced elements are 
required for resolution. In the remaining area between 
the body and bow shock, flow gradients are 
significantly smaller hence element spacing can oe 
larger. Since the flow structure is known apnon the 
flow field mesh is constructed to space elements 
closely together in the shock region ana boundary 
layer. The elements are also densely packeo along 
the stagnation line (y*0). In this region close to tne 
cylinder surface, the flow velocities are very small ana 
become zero at the stagnation point. Away from tne 
flow symmetry line, the subsonic flow accelerates 
nonuniformally to sonic speeds. 

The finite element model consists of 12,000 
quadrilateral elements in the flow domain ana 3,000 
quadrilateral elements for the cylinder with the same 
discretization along the flow-cylinder interface. A 
graded radial spacing was used in the cylinder. About 
35% of the fluid elements lie in the 0.02 inch thick 
boundary layer. The boundary layer mesh is 
graduated normal to the cylinder surface by an 
incremental factor of 1.1. The smallest element is at 
the flow stagnation point and is only 0.0004 inch long 
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in the direction normal to the cylinder surface. This 
extremely small element dimension is necessary to 
capture the temperature gradient and hence the 
aerodynamic heating rate. 

Typical flow solutions are shown in Figs. 4-7. 
Predicted density contours are compared in Fig. 4 with 
an experimental schlieren photograph 6 which is also a 
measure of flow density. The comparison of the shock 
shape and position indicates the global flow field is 
reasonably well predicted. Fig. 5 shows the excellent 
agreement between the computed aerodynamic 
pressure distribution along the cylinder surface and 
the experimental data 6 . The predicted and 
experimental pressures are normalized by their 
respective stagnation point pressures which were 
within 5% of each other. 

The major difficulty in the fluid analysis is the 
prediction of the aerodynamic heating rates because a 
very accurate resolution of the flow temperature 
gradient normal to the cylinder surface is required. 
The computed flow temperature distribution (Fig. 6) 
along the flow symmetry line (y=0) clearly illustrates 
the sharp gradients that must be resolved. The free- 
stream temperature increases abruptly from 435°R to 
about 3,900°R across the bow shock. Within a very 
thin layer at the flow stagnation point, the temperature 
drops sharply from 3,900°R to the surface temperature 
of 530°R. This steep temperature gradient within the 
thin layer, which is about 3% of the shock layer 
thickness, produces a high stagnation point heating 
rate. The predicted bow shock location compares well 
with the location predicted by empirical formulae 12 
and the schlieren photograph (Fig. 4). The predicted 
heating rate distribution normalized to the stagnation 
point heating rate at time zero compares well with the 
experimental results 6 as shown in Fig. 7. The 
predicted stagnation point heating rate (42.5 Btu/ft 2 - 
sec) is in excellent agreement with the Fay and Riddell 
solution 13 (42.5 Btu/ft 2 -sec) and a viscous shock layer 
solution 14 (41 .4 Btu/ft 2 -sec) but is quite lower than the 
experiment 6 (61.7 Btu/ft 2 -sec). The difference 
between the predicted and the experimental heating 
rate has not been resolved, but it could result from high 
temperature effects not taken into account by the 
analysis or from free stream turbulence in the test 
stream. 

Since the flow field reaches equilibrium 
instantaneously compared to the thermal response of 
the cylinder, this aerodynamic analysis was updated at 
two second intervals. The cylinder temperature 
contours at two seconds are shown in Fig. 8. The 
figure shows the maximum temperature of about 
700°R occurs at the stagnation point whereas the 
temperature along the back side of the cylinder 
remains at the ambient temperature of 530°R. The 
structure responds instantly to the pressure and 
temperature field hence a quasi-static structural 
analysis is adequate to predict the cylinder 
deformations and thermal stresses. The structural 
boundary conditions and the resulting circumferential 
thermal stress distribution on the cylinder are shown in 
Fig. 9. A maximum compressive stress of 32 ksi occurs 


on the cylinder outer surface near the stagnation point 
due to the high temperature gradients in that region 
The cylinder deformation is negligible compared to me 
cylinder diameter due to the relatively low temperature 
change and hence has a neglible effect on trie 
aerodynamic flow. However, the change in the surface 
temperature has a moderate effect on the flow fieio 
The aerodynamic heating rate is reduced in the region 
where the surface temperature has increases os 
shown in Fig. 7, due to a thickening of the thermal 
boundary layer and a lower fluid temperature gradiem 
at the cylinder surface. The stagnation point heating 
rate decreased nearly 8% from the initial level witn the 
cylinder temperature at 530°R. 

The coupling in this problem is moderate to wea* 
and occurs primarily between the flow ana thermal 
analyses. Since the flow and stress fields are no* 
coupled and develop instantaneously compared to me 
structural temperature field, an iterative solution 
between disciplinary solutions would be adequate O* 
course, the integrated analysis significantly easea aata 
transfer between disciplinary analyses. 

A major difficulty encountered in the flow fieia 
solution was the convergence of the aeroaynamic 
heating rate. The aerodynamic heating rate 
especially in the flow stagnation zone, convergea 
slowly because the small mesh needed to capture m« 
high temperature gradient required a very small time 
step to satisfy the computational stability criterion 
Hence, coupled with low flow velocities, perturbations 
take a large number of time steps to convect out o« the 
flow domain. 

To increase the convergence rate of tne 
aerodynamic heating and reduce the analysis 
computational time, the current analysis approacn 
must be modified. This may be done by applying an 
implicit solution technique for the flow oounaary iayet 
region to alleviate the computational stability constraint 
and thus increase the convergence rate Tme 
computational intensity can be further reduceo ov 
reducing the number of unknowns (i.e. gria points ana 
elements). This is achieved by adapting tne 
computational domain to the physics ot me 
phenomena creating an unstructured computational 
domain. 

These modifications were developea recently at 
the Aerothermal Loads Branch, NASA Langley 
Research Center, in a program called LARCNESS 16 
(Langley Adaptive flemeshing £ 0C j e an a iiavi£r 
S_tokes S.olver). The approach employs an 
implicit/explicit upwind finite element algorithm 
combined with an adaptive unstructured mesn 
refinement technique. The approach has been 
applied with excellent results to the flow problems 
described herein 15 . To demonstrate an application o> 
the adaptive unstructured mesh refinement technique 
a more complex flow behavior representing the 
second example where an oblique shock intersects 
the cylinder bow shock is described in Fig. 10(a). The 
intersection creates a supersonic jet that impinges on 
the cylinder 6 . The adaptive unstructured finite element 
model is shown in Fig. 10(b). The elements are 
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adaptively clustered to accurately capture the shock 
interaction phenomena and the flow solution in the 
boundary layer. The adaptive procedure, based on 
density gradients in the flow field, significantly reduced 
the number of grid points. This more complex flow 
field was resolved with 9000 grid points compared to 
12000 grid points used in this study for the simpler 
symmetric flow field. Similar adaptive procedures 
could be developed for the thermal and structural 
analysis. However, since the thermal analysis is 
transient the mesh would have to move with time to 
track the thermal response. The aerodynamic heating 
rate along the cylinder surface is compared with 
experimental results in Fig. 10(c). This new analysis 
approach for predicting complex flow behavior will be 
implemented as an alternative fluid module in LIFTS. 

CfloyectlyfilyXoalfid Leading Edqa 

A 0.25 inch diameter leading edge is used in the 
second example to study and identify the issues 
related to the thermal-structural analysis and the 
design of convectively cooled leading edges subjected 
to intense aerodynamic heating. The example 
represents the acceleration of a hypersonic vehicle 
through Mach 16 which causes the vehicle nose bow 
shock to sweep across the engine cowl leading edge 
from an outboard to an inboard position. The 
sequence of events is classified into the three 
conditions shown in Fig. 11. For condition I, the 
oblique shock, representing the vehicle nose bow 
shock and created by the Mach 16 free stream flow, 
passes outboard of the cowl leading edge. At this 
condition, the cowl leading edge is exposed to the 
Mach 3 flow behind the Mach 16 shock wave. It 
should be noted that the aerodynamic heating rate on 
the leading edge at this Mach 8 flow (condition I) is 
higher than that at the Mach 16 flow (condition III) due 
to the additional compression of the flow as it passes 
through the oblique shock wave. As the vehicle 
continues to accelerate through Mach 16, the vehicle 
oblique shock wave moves across the cowl leading 
edge (condition II) intersecting with the leading edge 
bow shock to produce transient shock wave 
interference heating. This interference heating results 
in a significant amplification of the heating rate on the 
leading edge as shown at typical times in Fig. 12. The 
envelope for the peak values and the distributions are 
idealizations of the experimental distributions given in 
Ref. 16. This idealization was used because the 
Taylor-Galerkin explicit flow algorithm would have 
been prohibitively expensive to perform the analyses. 
The interference heating rate reaches a peak value 
about 22° below the horizontal centerline of the 
leading edge. The oblique shock is assumed to move 
at a speed of two inches per second and passes the 
leading edge in 0.125 seconds. After the oblique 
shock passes, the leading edge is heated by the Mach 
16 flow (condition ill). 

The leading edge geometry, boundary conditions, 
and a schematic of the finite element thermal-structural 
model are shown in Fig. 13. The leading edge outer 
surface is subjected to the transient aerodynamic 
heating given in Fig. 12 and emits radiant energy to 
space. The inner surface is convectively cooled by the 


direct impingement of a sonic hydrogen jet stream with 
an inlet temperature of 50°R. The hydrogen flow rate 
is 0.104 Ibm/sec and the coolant film coefficient is 8 
BTU/ft 2 s. The finite element model consists of 1440 
quadrilateral elements for the leading edge and 90 
mass-transport conduction/surface convection 
elements representing the hydrogen coolant flow. The 
leading edge mesh is graded in the radial direction 
similar to the cylinder mesh shown in Fig. 3 but is 
uniform in the circumferential direction. An adaptive 
unstructured mesh based on the temperature 
gradients would significantly reduce the number of grid 
points and elements. 

Nickel-200 was first selected as the leading edge 
material because of its high ductility, high melting 
temperature, relatively low thermal expansion, and 
ease of fabrication. The nickel material properties 
used in the analysis are temperature dependent and 
can be found in Ref. 1 7. Since the structure will be 
brought to thermal equilibrium at condition I, the initial 
structural temperature was assumed to be equal to the 
coolant temperature of 50°R. Figure 14 shows the 
structure and coolant temperature histories at the 
selected locations shown in the inset during the 
sequence of events. At the early time, the leading 
edge is first subjected to the Mach 8 flow (condition I) 
for 0 s < t< 0.06 s. The figure shows that the steady- 
state condition is reached at the end of this time 
interval and the maximum temperature of about 700°R 
occurs at the horizontal centerline on the outer surface 
of the leading edge (location A). The inner surface 
(location B) remains at a relatively low temperature of 
350°R due to heat convection to the hydrogen coolant. 
As the hydrogen coolant flows along the coolant 
passage, its temperature increases from the inlet 
temperature of 50°R to 77°R at the exit (location C). 

During the transient interference heating (condition 
II, 0.06 s < t < 0.185 s), the outer surface temperature 
at the centerline (location A) peaks at 1500°R at the 
time of the oblique shock impingement (t = 0.12 s). 
The maximum surface temperature of approximately 
2500°R occurs at 22° below the centerline (location D) 
due to the maximum localized heating at i»0.14s. 
The predicted peak temperature is close to the melting 
temperature of the nickel indicating the current design 
concept is not acceptable. As the oblique shock 
passes ( t > 0.185 s), the leading edge is heated by the 
Mach 16 flow (condition III). The leading edge 
temperature drops due to the lower aerodynamic 
heating at this condition. 

To reduce leading edge temperatures which occur 
during the sequence of events just described, internal 
fins can be added to the convectively cooled area to 
increase the amount of energy absorbed by the 
coolant. Materials with "super* thermal conductivity 
can be used to reduce thermal gradients through the 
thickness and increase circumferential diffusion of the 
heat load. Copper and beryllium are the candidate 
materials because they provide higher thermal 
conductivity than nickel as shown in Fig. 15. The 
thermal conductivity of these two materials are 
reported to be exceptionally high at low temperature 18 . 
The peaks of the thermal conductivity for copper and 
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beryllium, which occur at 20°R and 70°R, are 
approximately 300 and 80 times higher than nickel, 
respectively. However the increased conductivity 
occurs over a very small temperature range offering a 
significant challenge to the computational procedure. 

The analysis is first repeated for a beryllium 
leading edge to investigate the reduction of the peak 
temperature response using a super thermal 
conductivity material. Internal tapered fins (-10 fins 
per inch) around the circumference of the leading 
edge are used to increase the convectively cooled 
area by a factor of five. However, an elementary 
analysis for the fin heat transfer efficiency shows that 
the fins increase the effective convection area by a 
factor of 2.5. In the analysis, the fins are taken into 
account by increasing the surface convection 
perimeter (P in eq. (5)) by a factor of 2.5. The 
circumferential conduction in the fin and the fin 
structural stiffness are not taken into account. 

The beryllium leading edge temperature history for 
the same locations previously shown for the nickel 
leading edge are shown in Fig. 16. The temperature 
contours at the Mach 8 flow condition I and at the time 
of the maximum interference heating are shown in Fig. 
17. Fig. 17(a) indicates that the entire leading edge 
acts as a super thermal conductor (illustrated by 
shaded area) with very low temperature rise at the 
Mach 8 flow condition I. As the oblique shock moves 
across the leading edge, the interference heating 
intensity increases, particularly in the shock 
impingement region, resulting in the localized high 
temperature zone shown in Fig. 17(b). This high 
temperature locally degrades the thermal conductivity 
and thus reduces the capability of conducting the high 
local heat flux from the leading edge outer surface to 
the hydrogen coolant. A maximum temperature of 
about 1000 °R occurs at the same time and location as 
for the case of the nickel leading edge. Details of the 
temperature response throughout the sequence of 
events indicate that the entire beryllium leading edge 
acts as a super thermal conductor at all times except 
during the short period of very intense aerodynamic 
heating shown in Fig. 17(b). No difficulties were 
experienced with the analysis procedure due to the 
extreme thermal property variation. 

Since copper has a higher thermal conductivity 
than beryllium at temperatures above 300°R, the 
analysis for the copper leading edge is performed to 
investigate further reduction of the maximum 
temperature at the time of the peak aerodynamic 
heating. The temperature history for the same 
locations previously shown for the nickel and beryllium 
leading edges are shown in Fig. 18. The maximum 
temperature is reduced to 766°R. Details of the copper 
leading edge temperature distributions as the shock 
sweeps across the leading edge are shown at the 
selected times in Figs. 19(a)-(d). it is interesting to 
note that under the Mach 8 flow condition, the 
maximum temperature of the copper leading edge 
(207°R in Fig. 19(a)) is higher than the beryllium 
leading edge (183°R in Fig. 17(a)). This is because 
copper has a lower thermal conductivity than beryllium 
in the narrow range of temperature near 200°R. 


The most severe thermal load on the copper 
leading edge, which occurs at the time of the peak 
aerodynamic heating (t=.l4 s in Fig. 19(c)), is used in 
the structural analysis to study the leading edge 
response. Since the thermal and structural models 
were indentical the transfer of nodal temperatures to 
the structural analysis was straight forward. 
Temperature dependent material strength properties 
for copper 17 are included in the analysis. A hydrogen 
coolant pressure of 1000 psi is assumed to apply 
uniformly along the leading edge inner surface. Along 
the outer surface, the aerodynamic pressure 
distribution is assumed to be the same shape as the 
aerodynamic heating 16 with a peak pressure of 1000 
psi. These pressure distributions and the structural 
boundary conditions are superimposed on the 
predicted leading edge structural response shown in 
Fig. 20. The figure shows the circumferential stress 
distribution on the deformed leading edge. The 
deformed leading edge shown in this figure is greatly 
exaggerated to highlight the detail of the deformed 
shape. At the oblique shock impingement location 
where the peak temperature and temperature 
gradients occur, the deformation is 0.0005 inches and 
the maximum compressive stress is 25 ksi. In other 
regions of the leading edge, the circumferential stress 
is moderate varying from 10 ksi in tension to 10 ksi in 
compression. 

The structural analysis was performed in two 
dimensions using the plane strain assumption. 
Although this is a conservative assumption, it clearly 
indicates the extreme axial compressive stresses 
(-100 ksi) that occur in the leading edge due to the 
extreme circumferential temperature gradient. The 
analysis suggests that the axial stresses for a three 
dimensional model of the leading edge may exceed 
the elastic limit of the material indicating a requirement 
for a nonlinear structural analysis. The current 
algorithm, which was derived for the inherently 
nonlinear fluid analysis, may be extended to three- 
dimensional nonlinear structural behavior with 
appropriate constitutive relationships. The ability to 
handle temperature dependent material properties 
and large deformations for two-dimensional problems 
has already been demonstrated 19 . The explicit 
algorithm with its restrictive stability criterion increases 
computational effort for quasi-static structural analysis, 
hence development of an acceleration procedure or a 
more efficient time marching scheme is needed. 

The structural response shows that the peak 
circumferential stress occurs at the oblique shock 
impingement location where the temperature and its 
gradients are maximum. The very high temperature 
and stresses could result in the failure of the leading 
edge, hence a more sophisticated structural analysis 
may be needed to include the capability to predict the 
permanent localized deformation including time 
dependency effects, such as cyclic loading and creep, 
as the oblique shock moves across the leading edge. 
LIFTS currently has the capability for the loading part 
of the cycle but not the unloading. The integrated 
approach of the LIFTS analyzer reduced the time and 
effort to perform the complete fluid-thermal-structural 
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analysis of the leading edges. Of course, in cases of 
strong interdisciplinary coupling the reduction in the 
effort would be even more substantial. 

Concluding Remarks 

A finite element approach for integrated fluid- 
thermal-structural analysis of aerodynamically heated 
leading edges was presented. The Navier-Stokes 
equations for high speed compressible flow and the 
associated thermal-structural equations were solved 
using a single finite element approach in one 
integrated, vectorized computer program called LIFTS 
(Langley Integrated £luid Ihermal Structural analyzer). 

The approach was used to study and validate the 
prediction of the flow behavior and the aerothermal 
loads by comparing with experimental data for Mach 
6.47 flow over a three inch diameter stainless steel 
cylinder. Comparison with experimental pressure and 
heat flux distributions is excellent. The thermal and 
structural analyses for the cylinder were performed 
and the response was examined to investigate 
benefits of the integrated interdisciplinary analysis and 
coupling. The coupling caused by the cylinder 
deformation that could alter the flow field was 
negligible. The fluid-thermal coupling from the 
increase in the cylinder temperature, which results in a 
reduction of the aerodynamic heating, was moderate 
but significant. The procedure automatically applies 
the true interface conditions (boundary conditions in a 
non-integrated procedure) as the solution is marched 
in time, thus avoiding the time and effort to interpolate 
data between different disciplinary models. The 
aerodynamic heating using the explicit Taylor-Galerkin 
algorithm converged very slowly because of the small 
time steps required for computational stability. 
Adaptive unstructured remeshing is required to reduce 
the computational effort. A new approach with a more 
efficient aerodynamic fluid algorithm and adaptive 
unstructured remeshing was described and will be 
incorporated into LIFTS as an alternative approach for 
the fluid analysis. 

Critical issues for the thermal-structural analysis 
were demonstrated by the analyses of a hydrogen 
cooled leading edge subjected to intense 
aerodynamic heating. The use of super thermal 
conductivity materials alleviates the elevated leading 
edge temperatures and gradients and hence the stress 
levels. The approach clearly demonstrates the ability 
to perform transient nonlinear thermal analysis for 
problems with highly nonlinear temperature 
dependent material properties and various boundary 
conditions. Under severe aerothermal loads which 
can cause intense localized leading edge 
temperatures, the present analysis indicates the need 
for more efficient algorithm with capability to predict the 
nonlinear, time dependent structural response. 

The applications have demonstrated the 
fundamental capability of the integrated analysis 
approach and provided insight into the features of the 
hypersonic flow over the leading edge, the leading 
edge thermal-structural response, and their 
interactions. Future analyses will study the flow- 


thermal-structural interactions for a more complex flow 
field that is altered by a larger structural deformation to 
further demonstrate the applicability and benefits of the 
approach for coupled interdisciplinary problems. 

Refercncfia 

1 . Shideler, J. L, Webb, G. L. and Pittman, C. M.: 
'Verification Tests of Durable Thermal Protection 
System Concepts,' Journal of Spacecraft and 
Rockets. Vol. 22, No. 6, November 1 985, pp. 598- 
604. 

2. Olsen, G. C., and Smith, R. E.: 'Analysis of 
Aerothermal Loads on Spherical Dome 
Protuberances," AIAA Journal. Vol. 23, No. 5, 

May 1985, pp. 650-656. 

3. Glass, C. E. and Hunt, L. R.: "Aerothermal Tests of 
Spherical Dome Protuberance on a Flat Plate at a 
Mach Number of 6.5," NASA Technical Paper 
2631, December 1986. 

4. Thornton, E. A. and Dechaumphai, P.: "Coupled 
Flow, Thermal and Structural Analysis of 
Aerodynamically Heated Panels," Presented at 
the AIAA/ASME/AHS 28th Structures, Structural 
Dynamics and Materials Conference, Monterey, 
California, April 6-8, 1987, AIAA Paper No. 87- 
0701. 

5. Thornton, E. A. and Dechaumphai, P.: "Finite 
Element Prediction of Aerothermal-Structural 
Interaction of Aerodynamically Heated Panels," 
Presented at the AIAA 22nd Thermophysics 
Conference, Honolulu, Hawaii, June 8-10, 1987 
AIAA Paper No. 87-1610. 

6. Wieting, A. R.: "Experimental Study of Shock 
Wave Interference Heating on a Cylindrical 
Leading Edge," Ph.D. Dissertation, Old Dominion 
University, Norfolk, Virginia, 1987, also NASA 
TM- 100484, 1987. 

7. White, F. M.: Viscous Fluid Flow . McGraw-Hill, 
1974. 

8. Thornton, E. A. and Wieting, A. R.: "Finite Element 
Methodology for Transient Conduction/Forced 
Convection Thermal Analysis," Progress in 
Astronautics and Aeronautic: Heat Transfer- 
Thermal Cont rol and Heat Pipes. Vol. 70, Edited 
by Walter B. Olstad, AIAA, New York, pp. 77-103. 

9. Bey, K. S., Thornton, E. A., Dechaumphai, P. and 
Ramakrishnan, R.: "A New Finite Element 
Approach for Prediction of Aerothermal Loads - 
Progress in Inviscid Flow Computations," 
Presented at the AIAA 7th Computational Fluid 
Dynamics Conference, Cincinnati, Ohio, July 15- 
17, 1985, AIAA Paper No. 85-1533, also NASA 
TM 86434. 

10. Thornton, E. A. and Dechaumphai. P.: "A Taylor- 
Galerkin Finite Element Algorithm for Transient 
Nonlinear Thermal-Structural Analysis," 


9 


Presented at the AIAA/ASME/AHS 27th 
Structures, Structural Dynamics and Materials 
Conference, San Antonio, Texas, May 19-21, 
1986, AIAA Paper No. 86-091 1 . 

1 1 . Thornton, E. A., Dechaumphai, P. and 
Vermaganti, G.: "A Finite Element Approach for 
Prediction of Aerothermal Loads," Presented at 
the AIAA/ASME 4th Joint Fluid Mechanics, 

Plasma Dynamics and Laser Conference, Atlanta, 
Georgia, May 12-14, 1986, AIAA Paper No. 86- 
1050. 

12. Billig, F. S.: "Shock-Wave Shapes Around 
Spherical and Cylindrical Nosed Bodies," 

Journal of Sp acecraft and Rockets. Vol. 4, 

No. 6, June 1967, pp. 822-823. 

1 3. Fay, J. A. and Riddell, F. R.: "Theory of 
Stagnation Point Heat Transfer in 
Dissociated Air," Journal of the Aeronautic 
Sciences . Vol. 25, No. 2, February 1958, 
pp. 73-85. 

14. Holcomb, J. E., Curtis, J. T., Shope, F. L.: "A 
New Version of the CVEQ Hemisphere Viscous 
Shock Layer Program for Equilibrium Air," AEDC- 
TMR-85-V7, February 1985. 

15. Thareja. R. R., Stewart, J. R., Hassan, O., Morgan, 
K., and Peraire. J.: "A Point Implicit Unstructured 
Grid Solver for the Euler and Navier-Stokes 
Equations,” Presented at the AIAA 26th 
Aerospace Sciences Meeting, Reno, Nevada, 
January 11-14, 1988, AIAA Paper No. 88-0036. 

16. Holden, M. S., and Wieting, A. R., Moselle, J. R., 
and Glass, C.: "Studies of Aerothermal Loads 
Generated in Regions of Shock/Shock Interaction 
in Hypersonic Flow," Presented at the AIAA 26th 
Aerospace Sciences Meeting, Reno, Nevada, 
January 11-14, 1988, AIAA Paper No. 88-0477. 

1 7. Metals Handbook Committee: Metals Handbook . 
Eighth Edition, Americal Society for Metals, Ohio, 
1975. 

18. Touloukian, Y. S., Powell, R. W., Ho, C. Y., and 
Klemens, P. G.: "Thermal Conductivity for Metallic 
Elements and Alloys," Thermophysical Properties 
of Matter, Vol. 1, IFI/Plenum, New York, 1970. 

19. Dechaumphai, P., Wieting, A. R., and Thornton, E. 
A.: "Thermal Structural Performance of an 
Actively Cooled Leading Edge Subjected to Type 
IV Shock Wave Interference Heating," Third 
National Aero-Space Plane Symposium June 2- 
4, 1987, Paper No. 24. 



Fig. 1 Fluid-thermal-structural interactions on an 
aerospace plane scramjet engine leading 
edge. 



Fig. 2 Experimental configuration for flow over a 

three inch diameter cylinder in the 8-Foot High 
Temperature Tunnel. 
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Fig. 3 Fluid-thermal-structural finite element model 
for flow over a three inch diameter cylinder. 





Fig. 4 Comparison of density contours with 

Schlieren photograph for Mach 6.47 flow over 
a three inch diameter cylinder. 
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Fig. 5 Comparative surface pressure distributions 
on a three inch diameter cylinder subjected 
to Mach 6.47 flow. 
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Fig. 7 Comparative surface heating rate 

distributions on a three inch diameter cylinder 
subjected to Mach 6.47 flow. 



Fig. 8 Temperature contours on a three inch 
diameter cylinder at 2 seconds. 
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Fig. 6 Fluid temperature distribution along the 

centerline of a three inch diameter cylinder for 
Mach 6.47 flow. 



8 ksi 


Fig. 9 Circumferential stress contours on a three 
inch diameter cylinder at 2 seconds. 






a) Oblique shock and bow shock interaction 
on cylinder. 



b) Adaptive mesh with enlarged section. 
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c) Comparative surlace heating rate 
distribution. 
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Fig. 1 1 A 0.25 inch diameter hydrogen cooled leading 
edge with shock wave interference heating. 



Fig. 12 Leading edge transient interference heating 
due to vehicle acceleration at Mach 16. 



Fig 10 Adaptive mesh refinement flow model and 
aerodynamic heating for cylinder with 
impinging oblique shock and bow shock 
interaction. 


Fig. 13 A schematic thermal-structural finite element 
model of 0.25 inch diameter leading edge with 
boundary conditions. 
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Fig. 14 Nickel leading edge temperature response. 


a) Mach 8 flow. 
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b) Maximum interference heating. 


Fig. 15 Temperature dependent thermal conductivity Pig. 17 Temperature contours and super thermal 

for copper, beryllium, and nickel. conductivity region for beryllium leading edge. 
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Fig. 1f> Beryllium leading edge temperature response. Fig. 18 Copper leading edge temperature response. 
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